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Results from numerical simulations of full, two flavor QCD thermodynamics at N t = 4 with 
domain wall fermions are presented. For the first time a numerical simulation of the full QCD phase 
transition displays a low temperature phase with spontaneous chiral symmetry breaking but intact 
flavor symmetry and a high temperature phase with the full SU(2) x SU(2) chiral flavor symmetry. 
It is found that close to the transition in the high temperature phase the U(1)a axial symmetry 
is broken only by a small amount. This result is of particular interest because of the connection 
between U{\)a symmetry breaking and the expected order of the transition. 
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On ■ 

' The upcoming RHIC experiments will reproduce the extreme conditions that existed in the early universe and it is 
hoped that the predicted QCD phase transition from hadronic matter to a quark-gluon plasma will be observed. The 
i non-perturbative nature of these phenomena makes precise theoretical predictions difficult. Lattice gauge theory pro- 
vides a first-principles approach to QCD thermodynamics and, in principle, a way to perform precise non-perturbative 
calculations of thermodynamic quantities at equilibrium. Quantities such as the value of the critical temperature, the 

£N| ■ width of the critical region and the order of the transition can then be used by phenomenological models to make 
contact with experiment. While great progress has been made since the inception of the lattice regulator 25 years 

^vq , ago there remain very significant uncertainties that prevent these calculations from being viewed as unambiguous 
predictions of QCD. Since these calculations are computer simulations one could think that progress in this field is 
tied with progress in computer speeds. While this is partially true, it is not the real culprit especially given the fast 

Q\ . growth in computer speeds during this period. The true reason lies at the most simple and fundamental level of the 

O^ \ lattice theory. 

When the first derivative of the fermion kinetic energy term is made into a lattice difference operator the number 
'"" J" 1 ■ of light species multiplies by a factor of 2 d where d is the dimension of space-time. This is the well known doubling 
^H, problem and it turns out that at non-zero lattice spacing in order to remove the extra species and still have a local 
Lagrangian one has to compromise the global flavor chiral symmetries of the theory |2J . There are two popular methods 
to put fermions on the lattice: Wilson and staggered fermions. Both of these methods at finite lattice spacing a break 
the chiral symmetry. The chiral symmetry is recovered together with the Lorentz symmetry as the continuum limit 
a — > is taken. And here lies the problem. In a numerical simulation a decrease of the lattice spacing by a factor 
of 2 requires an increase in computation by a factor of 2 8-10 depending on the parameters. In order to control the 
amount of chiral symmetry breaking induced by the regulator prohibitively large computing resources are required. 
Clearly faster computers although one day may be able to overcome these large factors can not be the answer. 

A way out of this forbidding state of affairs came from D. Kaplan || through an unexpected avenue. He formulated 
a new fermion lattice regulator with the name domain wall fermions (DWF). The first numerical simulations using 
DWF in vector gauge theories were done in Q and were followed by §, @, 0, §, §, [0. In this work the latest 
results from simulations of dynamical QCD thermodynamics are presented. Preliminary results of this work have 
already appeared in || , |lC| . For promising alternatives to domain wall fermion simulations see Jtl| . 
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II. DOMAIN WALL FERMIONS. 



Domain wall fcrmions introduce an extra direction of space-time with size L s . The fermion fields are five-dimensional 
while the gauge fields remain four-dimensional and couple the same way to all extra fermion degrees of freedom. 
The boundary conditions along the fifth direction are free (domain wall) and although the extra fermion degrees of 
freedom are heavy a light Dirac fermion surface mode develops on the boundaries with its positive chiral components 
exponentially bound on one wall and its negative components on the other. In the L s — > oo limit the mode is a massless 
Dirac fermion. At finite L s the residual mixing introduces an exponentially small mass. Since these properties are 
maintained at any lattice spacing DWF provide a way of separating the recovery of the chiral symmetry from the 
recovery of the Lorentz symmetry. Now at a fixed lattice spacing the chiral symmetry can be recovered by increasing 
L s . Unlike Wilson or staggered fermions the computing cost in recovery of the chiral symmetry is only linear in L s \ 
These remarkable properties may bring within reach physically realistic studies of lattice QCD thermodynamics. 

The subject of DWF has a large volume of work behind it. The reader is referred to Q and references therein for 
an introduction. Here, the version of DWF as presented in jl2| with the modifications in Q is used. The theory 
has five parameters: the four-dimensional lattice volume V , the inverse gauge coupling squared (3, a bare mass that 
explicitly mixes the two chiralities bound on the domain walls to/, the number of sites of the fifth direction L Sl and a 
five-dimensional mass that can be thought of as the domain wall height too. The parameters To/, L s and too set the 
bare quark mass. In free theory the bare quark mass m e s is given by [Q: 

m c g = to (2 - too) [to/ + (I - m ) Ls ] , < to < 2 (1) 

There are many important theoretical issues regarding DWF. Perhaps the most crucial one relates to the localization 
of the chiral modes in the interacting theory. Does the localization remain exponential as in free theory or does it 
become power law? How does the decay rate depend on the lattice spacing? These questions were addressed in detail 
in Q in the context of the two flavor Schwinger model. There it was found that the decay is exponential and that 
the decay rate becomes faster as the continuum limit is approached. The last result can be understood by a simple 
examination of eq. [I]. In the interacting theory mo will be renormalized and therefore optimum localization will be 
shifted away from too = 1. However, even if too is adjusted to absorb this renormalization, optimal decay will not be 
possible since the optimal value for too will fluctuate from one background gauge field to the next. At strong coupling 
were the gauge field fluctuations are large the corresponding fluctuations in the optimal value of too will also be large 
and the decay rate will become slow. This scenario is reminiscent of the way the bare mass of Wilson type fermions 
is used to cancel the renormalization of the interacting theory. However, because the massless point changes from 
configuration to configuration masslesssness can only be achieved up to order of the lattice spacing. Here the same 
is true but the fluctuations in mass are raised to the L s power and therefore can be de-amplified by any amount for 
large enough L s . 

There is another related source of de-localizing effects. It has been shown that the transfer matrix along the extra 
direction for certain gauge field configurations can develop unit eigenvalues . For such configurations the surface 
states completely de-localize. Although, the set of such configurations is of zero measure, configurations in their 
vicinity will produce poor localization. These configurations are special to the lattice since they allow the passage 
from one topological sector to another. Again, as the continuum limit is approached these configurations will have a 
negligible Boltzmann weight and their de-localizing effects will diminish. 

Another issue relates to the positivity of the transfer matrix along the fifth direction. Details regarding this can be 
found in [T^ ]. Here, it suffices to mention that if L s is kept even and if all Green's functions are extracted from sites 
along the fifth direction that have a distance from the boundaries of an even number of sites the relevant transfer 
matrix is the square of the single step one and is therefore positive. 

Finally the number of light flavors of the free theory depends on too . For example if < too < 2 the theory has 
one light flavor but if 2 < too < 4 the theory has four light flavors. These too ranges are modified in the interacting 
theory and one should make sure that the needed number of flavors is produced M, [pi . 
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III. TWO FLAVOR SCHWINGER MODEL. 



The two flavor Schwinger model was simulated numerically in [J4J . If the chiral symmetries of the massless model 
are intact, the chiral condensate will have a zero vacuum expectation value, < WW >= at all couplings. Therefore, 
the deviation of < WW > from zero can serve as a measure of the amount of chiral symmetry breaking induced by 
the regulator. In figure [l] the results from these simulations are presented. 
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FIG. 1. The chiral condensate normalized by the photon mass for m/ = 0, mo = 0.9. From top to bottom the lines 
correspond to lattice spacing a of 1/6, 1/8, 1/10 and 1/12. The solid lines are fits to coe clLs . 

There are three observations that can be made from this figure: 

1) For each lattice spacing there are two different decay rates. In Q the fast rate was related to the physics of the 
zero topological sector while the slow rate was related to the physics of topology changing configurations. 

2) Fits to an exponential decay form c$e~ ClLs have x 2 per degree of freedom around one indicating that the localization 
of the light states in the interacting theory is consistent with exponential decay. Also, power law behavior can be 
excluded with some confidence for the a = 1/12 data since a power law fit has a % 2 per degree of freedom of 32. 

3) The decay rates become faster as the lattice spacing is decreased. This is a key feature that makes DWF a viable 
lattice fermion regulator. 

IV. QUENCHED QCD. 

Before DWF can be used in QCD an analysis similar to the one for the Scwhwinger model must be done in order to 
map the parameter space. Below the results from such an analysis using quenched QCD are presented. Preliminary 
results can be found in 0, ||. 

As mentioned in section || too controls the number of light flavors. Since this parameter is renormalized, it is 
important to understand which ranges correspond to different flavor sectors. In figure [2^, < WW > is plotted as a 
function of too- Three regions can be distinguished: < WW > is zero for too < 1, it has a value around 0.001 for 
1 < mo < 2.5 while it has a value that is four times larger for 2.5 < mo < 4.0 . Therefore, these three regions 
correspond to a zero, two and eight flavor theory (for algorithmic reasons only an even number of flavors is simulated; 
had it not been for this restriction the three regions would correspond to zero, one and four flavors). 

In order to investigate the L s dependence of the theory the pion mass m^ is measured for several values of to./ at 
fixed L s and mo for (3 — 5.7 on an 8 3 x 32 lattice. Its value is then extrapolated to zero quark mass. This exercise 
is repeated for several values of L s . The mj = extrapolated value as a function of L s is shown if figure |b for 
mo = 1.65. The solid line is a fit to Co + c\e~ C2Ls . The fit has x 2 /d.o.f. « 1 indicating that the localization is 
consistent with exponential decay. The constant part of the fit is due to finite volume. On a finite box of size L and 
periodic boundary conditions the longest correlation length is L/2 and therefore the smallest mass should be ss 2/L. 
In this case m w (mf = 0, L s = oo) w 0.22 m 2/8. That this is indeed a finite volume effect was supported by a second 
simulation on a 16 3 x 32 lattice. 
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FIG. 2. (a) < tyty > extrapolated to my ■ — vs. mo for various L s . (b) extrapolated to m/ = vs. L s . 

Finally, the rho and nucleon masses were measured for several to/ and L s . Their to/ — extrapolated values 
showed almost no dependence on L s for 10 < L s . The nucleon to rho mass ratio was found to be jm p — 1.508(65) 
to be compared with the physical value 1.221(2). The discrepancy is expected since the lattice spacing is rather 
large, a w 0.2/to, and the box size only marginally large 1.6/to. Given the modest volume size these results compare 
favorably with previous quenched QCD results using Wilson and staggered fermions at larger volumes. 

V. DYNAMICAL QCD THERMODYNAMICS. 

With this ground work in hand dynamical QCD thermodynamics can be investigated. The first exploratory work 
is done on small lattices of size 8 3 x 4. Preliminary results have been presented in ||. In figure ||a < tyty > as a 
function of (3 is shown for a fixed L s and m/ at too = 1.9. A rapid crossover is observed indicating the presence of 
a phase transition. To confirm that indeed there are two phases < '5^' > is plotted in figure ||b as a function of to/ 
at L s — 16 for (3 = 5.20 below the transition and (3 = 5.45 above. Clearly below the transition < 'J' 1 !' > extrapolates 
to a non zero value indicating spontaneous breaking of the SU(2) x SU(2) chiral symmetry down to SU(2) flavor. 
Because the Langrangian is explicitly symmetric under the SU (2) flavor there are three degenerate Goldstone pions. 
Above the transition < 4"!' > extrapolates to a non zero but very small value (which could be further reduced by 
increasing L s ) indicating that the full SU(2) x SU(2) chiral symmetry is essentially restored. This illustrates the 
unique properties of DWF. 

To further study the properties of the dynamical theory < ^"5 > is plotted as a function of L s for to/ = 0.02, 
mo = 1.9 both below (f3 — 5.20) and above (/3 = 5.45) the transition in figure ^a. The fits are of the form co + cie~ C2Ls . 
The constant term is the value of the condensate at the non-zero mass to/. The % 2 /d.o.f. is small indicating that 
the localization of the dynamical theory is consistent with exponential. As seen in the Schwinger model the decay 
rate becomes faster as the lattice spacing is made smaller (increasing 0). At f3 = 5.45 the L s — oo value of < "I"!' > 
is already approached to within a few percent at L s = 16 while at j3 — 5.20 a larger L s is needed. Although these 
results are encouraging it must be stressed out that at this large lattice spacing other observables may have different 
decay characteristics. 

Finally it is important to investigate the behavior of the phase transition as a function of too. In figure |]b the 
Wilson line is shown as a function of (3 for various mo. A similar plot has also been obtained for < \E f \E f > but because 
the wavefunction renormalization of the fermion field depends on too, it is hard to display in a single figure. From 
left to right the figure corresponds to too = 2.15, 1.90, 1.80, 1.65. The qualitative features of these curves are similar 
but the critical value of f3 is a strong function of too. This is not surprising since /3 is a bare parameter. One observes 
that as to is reduced the critical [3 approaches the quenched value (3 C = 5.7 indicating that the regime of zero flavors 
is approached. The role of toq in thermodynamics should be studied with some caution. For example, from free field 
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studies the range of small momenta accessible to the light flavor is restricted if mo is set close to its lower bound for 
one flavor physics. Similar effects can occur if mo is set close to its upper bound for one flavor physics. Since such a 
restriction may affect the thermodynamics one should set mo somewhere in the middle of the allowed range. 
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FIG. 3. < 99 > around the QCD finite temperature phase transition. 
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FIG. 4. (a) < 99 > vs. L a below and above the transition, (b) The Wilson line \W\ vs. /3 for various mo. 



VI. THE U(1) A SYMMETRY JUST ABOVE THE QCD DECONFINING TRANSITION. 



If the U(1)a symmetry above but close to the deconfining phase transition is broken the transition is expected to 
be second order while if it is not broken or only softly broken the transition may be first order |H| . This important 
issue can be investigated by direct lattice QCD simulations just above the transition. Using staggered fermions (l(J 
it was not possible to draw unequivocal conclusions mainly because staggered fermions may not produce the zero 
mode effects necessary for anomalous breaking of U(1)a 0], |@- On the other hand DWF are known to posses 
exact and robust zero modes in the L s — oo limit Jf3| . In |f9| ] it was shown that for classical instanton backgrounds 
these properties are maintained with high accuracy even at finite L s rj 10. In order to confirm the presence of zero 
modes in a situation with quantum fluctuations < \E r \& > versus m/ from a numerical simulation of quenched QCD 
just above the transition is plotted in figure ||a. Since in the quenched approximation the fermionic determinant is 
set to one, the zero mode effects are not suppressed and < 99 > diverges as 1/m/. 

Having established that DWF exhibit the desired zero mode effects the U{\)a symmetry is probed by numerical 
simulations of full QCD using DWF on a large lattice above the deconfining transition. The difference m$ — m n of 
the screening masses of the delta and pion particles is used as a measure of anomalous symmetry breaking. The 
screening masses are measured from the exponential fall-off of the relevant two point Green's function along a spatial 
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direction. Since the delta and pion Green's function are related by a U transformation ms — m n should be zero in 
the zero mass limit if the symmetry is not broken and non-zero otherwise. In figure ms — m^ versus m/ is shown 
for ft = 5.45 and ft = 5.40 on a 16 3 x 4 lattice with L s = 16 and mo — 1.9. The critical ft is around 5.325. The lines 
are fits to cq + c^m 2 and have x 2 /d.o.f . « 1. The absence of a linear term indicates that chiral symmetry is effectively 
restored. The mj = extrapolated values are 0.087(17) at ft = 5.40 and 0.031(9) at ft = 5.45. Although both are 
not zero their value is small when compared with ms and m, which are « 1.3. Universality arguments require that if 
the QCD phase transition is to be second order, the anomalous U(1)a must be broken. It is an open question as to 
whether the small size of the U(1)a symmetry breaking seen here is sufficient to support this theoretical prediction 
that the two- flavor QCD phase transition is second order |La|. 




FIG. 5. (a) < tyty > in quenched QCD just above the transition, (b) ms — m w in dynamical QCD just above the transition. 



VII. THE CHARACTER OF THE QCD PHASE TRANSITION. 

In order to determine the critical temperature and investigate the width of the critical region and the order of the 
transition, simulations of full QCD with DWF on 16 3 x 4 lattices close to the transition are performed for mo = 1.9, 
rrif = 0.02 and L s = 24. The chiral condensate is plotted as a function of ft in figure ||. For each (3 two separate 
simulations are done one with an ordered initial configuration (diamonds) and one with a disordered one (crosses). 
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FIG. 6. < > vs. 13 on a 16 3 x 4 lattice with m = 1.9, m f = 0.02 and L s = 24. 

A zero temperature dynamical simulation on an 8 3 x 32 lattice at ft = 5.325 with mo = 1.9 mf — 0.02 and L s = 24 is 
done in order to set the scale close to the transition. The valence extrapolated p mass is found to be m p — 1.198(57). 
If it is used to set the scale it results to a critical temperature T c « 161MeV. The pion mass at m/ = 0.02 is 
inn k 412MeV. Further analysis attributes this large value to finite L s effects. This is under current investigation. 

Finally, from figure ^ one can see that the ordered and disordered points at ft = 5.325 disagree even though each 
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was determined after 400 thermalization sweeps were discarded. These small discrepancies visible for j3 = 5.3 and 
5.325 suggest the presence of a critical region and indicate the need for collecting more statistics, a process currently 
underway. Also, must be reduced before the order of the transition can be identified. This is currently in progress. 



VIII. CONCLUSIONS. 



A novel lattice fermion regulator, domain wall fermions, was used to study QCD thermodynamics. For the first 
time, this regulator offers the possibility to control the notorious lattice chiral/flavor asymmetries that have impeded 
numerical studies of QCD thermodynamics. The U{1)a symmetry above the deconfming transition was found to be 
broken only by a small amount. The critical temperature for the transition was calculated and was found to be in 
agreement with previous estimates. Detailed studies of the transition region are currently underway with the hope 
that a consistent picture of the order of the transition and the relatively small U(\)a breaking will emerge. However, 
smaller pion masses are needed before any conclusions can be drawn and such work is now underway. 



IX. ACKNOWLEDGMENTS. 



The simulations were done on the 400 Gflops QCDSP supercomputer |2fj] at Columbia. This research was partially 
supported by DOE under grant # DE-FG02-92ER40699 and for PV partially under grant # NSF-PHY96-05199. 



K.G. Wilson, Phys. Rev D10 (1974) 2445. 

H.B. Nielsen and M. Ninomiya Nucl. Phys. bf B185 (1981) 20. 

D.B. Kaplan, Phys. Lett. B288 (1992) 342; Nucl. Phys. B30 (Proc. Suppl.) (1993) 597. 

P.M. Vranas, LATTICE 96, Nucl. Phys. B53 (Proc. Suppl.) (1997) 278; P.M. Vranas Phys. Rev. D57 (1998) 1415. 
T. Blum and A. Soni, Phys. Rev. D56 (1997) 174; Phys. Rev. Lett. 79 (1997) 3595. 

P. Chen, N. Christ, G. Flemin g, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
contribution to LATTICE 98, |hep-lat/980915£ . 

P. Chen, N. Christ, G. Flemin g, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
contribution to LATTICE 98, |hep-lat/981102e . 

P. Chen, N. Christ, G. Flemin g, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
contribution to LATTICE 98, |hep-lat/9811013 . 

P. Chen, N. Christ, G. Fleming, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
contribution to LATTICE 98, talk by A. Kaehler. 

P. Chen, N. Christ, G. Fle ming, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
contribution to ICHEP 98, |hep-lat/9812011 . 

H. Neuberger, Phys. Rev. D57 (1998) 5417; Phys. Lett. B417 1 41 (1998); Phys. Rev. Lett. 81 (1998) 4060; tiep-lat/990100^ ; 



[20] 



U.M. Heller, R. Edwards and R. Narayanan |hep-lat/9807017| ; [hep-lat/981103C| C. Lieu, [hep-lat/9811008 
V. Furman, Y. Shamir, Nucl. Phys. B439 (1995) 54. 
R. Narayanan, H. Neuberger, Phys. Lett. B302 (1993) 62; Phys. Rev. Lett. 71 (1993) 3251; Nucl. Phys. B412 (1994) 574; 
Nucl. Phys. B443 (1995) 305. 

P. Chen, N. Christ, G. Fleming, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
in preparation. 

R. Pizarski and F. Wilczek, Phys. Rev. D29 (1984) 338. 

S. Chandrasekharan, D. Chen, N.H. C hrist, W. Lee, R. Mawhinney, and P.M. Vranas, CU-TP-902, [hep-lat 9807018 . 
J.B. Kogut, J.F. Lagae, D. K. Sinclair, |hep-lat/9801020t 
A. Kaehler in preparation. 

P. Chen, N. Christ, G. Fleming, A. Kaehler, C. Malureanu, R. Mawhinney, G. Siegert, C. Sui, P. Vranas, and Y. Zhestkov, 
Phys. Rev. D59 (1999) 054508. 

D. Chen, P. Chen, N. Christ, R. Edwards, G. Fleming, A. Gara, S. Hansen, C. Jung, A. Kaehler, A. Kennedy, G. Kilcup, 
Y. Luo, C. Malureanu, R. Mawhinney, J. Parsons, C. Sui, P. Vranas, and Y. Zhestkov, contribution to LATTICE 98, 



hep-lat/9810004 



7 



